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THREE-DIMENSIONAL BOUNDARY LAYER ANALYSIS PROGRAM 
"BLAY" AND ITS APPLICATION 

*/207 


Ken-ichi Matsuno and Tomiko Ishiguro 


1 . Introduction 

Among the important elements in the transonic aerodynamic 
design of a wing is the evaluation of viscous strength. The flow 
field of the transonic region itself is an aerodynamically critical 
region. For example, even when considering only the shock wave on 
the wing surface, it is a well known fact that the shock position 
and shock strength differ significantly when comparing the cal- 
culated results for the case where non-viscous flow is assumed and 
where viscous effects are taken into account. The viscous effects 
for adhesive flow under flight cruise conditions can be evaluated 
with sufficient accuracy even by applying boundary layer equations. 
Therefore, the computational aerodynamic technique currently being 
utilized in aerodynamic design involves the method of compensating 
in the transonic totai potential flow analysis program by incorpor- 
ating the boundary layer effects. This method is more realistic 
than considering developmental factors such as computing costs and 
the computer hard level calculation method. Thus, a variety of test 
runs are currently being conducted. 

The BLAY (boundary layer calculation program) is the acronym 
for a program that analyzes rigorously the three-dimensional bound- 
ary layer of a wing having an arbitrary planar configuration. It 
is based on a research program developed in the 1979-1980 period 
which subsequently was further developed into a general purpose 

* 
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program with the objective of its conversion into a production 
code . 


The purpose of this paper is to present a general survey of 
the development status among several nations in their three-dimen- 
sional boundary layer analysis programs dealing with wing configura- 
tions; to discuss the characteristics of the BLAY program in a com- 
parative evaluation against those programs; and to point out that, 
among all presently existing computation programs, the BLAY is out- 
standing in many factors including computation efficiency. As 
examples of computations employing the BLAY, we will introduce the 
results on some wing configurations of recent interest, including 
the forward swept wing. 

2 . Three-dimensional boundary layer analysis program on a wing 

Since we are here particularly interested in solving rigorously 
the boundary layer equation(s), the so-called integration method, 
based on Karman's integration equation, will be omitted from con- 
sideration on the ground that it is an approximation analysis method. 
The three-dimensional boundary layer analysis program which addresses 
wing configurations is being developed energetically in particular 
by airframe manufacturers in the United States, and by the national 
aeronautical research laboratories in Europe, as can be deduced 
from the news on those subject areas. These were reported continously 
in the 1977-1979 period as shown in Table 1. This is exactly the 
time frame during which Jameson's Transonic Total Potential Flow 
Analysis Program (FLO 22) was published. In our country, too, R6D 
was initiated in 1978 at the Aero Tech Lab (Ko Gi Ken) and followed 
the course shown in Table 1, and produced the current program, BLAY. 

2.1 Coordinate system structure and boundary layer equations 

*V' — * ~ 

In the case of computational programs of this type, which are 
evaluated against computation criteria for a real (finite) wing, 
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the coordinate system to be used may constrain the applicable scope 
of the program or. the region within which the computation can be 
performed. Figure 1 compares the coordinate system employed in 
several computation programs and the BLAY system. The coordinate 
system used by Nash and Scruggs, and others [2,10] is a polar coor- 
dinate system (r, 0, y) applicable to the straight line portion of 
a tapered wing and cannot be applied to parts with strakes or kinks. 

In MacLean's [4] coordinate system, a cross flow coordinate system 
is created from the leading edge stagnation line and can thereby be 
applied to various wing configurations, but has the following defects: 
generally, there occurs a region on the wing tip(s) where computation 
is impossible; it is difficult to establish initial conditions at 
the wing root; the computation points at the wing root differ at all 
coordinate lines and processing becomes complex; the non-viscous data 
necessary as a boundary condition are generally given at fixed per 
cent chord,, and fixed per cent span. Therefore, there is no match 
with the coordinate points of the boundary layer computation, and it 
becomes necessary to employ a complicated interpolation routine; 


TABLE 1. Three-dimensional boundary layer analysis programs for 
wings (references 1 through 13) 

Boundary Layer Methods for Finite Wings 


United States 
Lockheed Corp. 

Nash and Scruggs 
Nash and Scruggs 

(1972) 

(1978) 

Explicit Scheme 
Implicit Scheme 

McDonnel 

Cebeci .Kaups and Ransey 

(1977) 

Regular Box Scheme 

Douglas Corp. 


(1978) 

Modification of a scheme 
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/Zig-zag Box Scheme 



(1979) 

VCharacteri stic Box Scheme 


Boeing Co. 

McLean 

McLean and Randall 

(1977) 

(1973) 

Implicit Scheme 
Pilot code 

Europe 
Bri tain 

Bradshaw, Mizner 

and Unswcrth 

(1975) 

for straigt tapered wings 

Netherlands 

Lindhout and De 
L indhout et al . 

Boer 

(1976) 

(1979) 

for laminar flows 
for turbulent flows 

West Germany 

Kordul 1 a 
Rostgi and Rodi 


(1977) 

(1979) 

Krause Ziy-zag Scheme 
for straight tapered wings 

Japan NAL 

(1978 : start ) 





(1979) Research code ; for rectangular wings 

Presentation of a new scheme ( PC-CN Scheme ) 
(1980) Research code ; for arbitrary finite wings 
(19S1) Modification of the scheme 
(193?) Development of. a pilot code 
(1983) BLAY code 
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similarly, in the case of output, a need arises to redistribute the 
data on a fixed per cent span. Among the coordinate systems used 
by Kordulla [9], the coordinates for the systems documented in the 
cartesian system are normalized in the chord direction at 0 at the 
leading edge and at 1 (unity) at the trailing edge, and is appro- 
priate for input/output of data at a fixed per cent chord and fixed 
per cent span. Furthermore, in this coordinate system, the velocity 
component is taken directly in the direction of the cartesian system. 


The coordinate system by Cebeci, et al. [8] uses non-cross flow 
coordinates that establishes fixed per cent chord ;and fixed per cent 
span as coordinates (x,y) initially. In this case, the number of 
terms appearing in the boundary layer equation is somewhat greater 
than that for the cross flow system; and is characterized by the 
fact that the establishment of initial conditions is direct, and, 
therefore, easy, at the leading edge stagnation line, at the wing 
root area and at the wing tip area. Also it is suitable for data 
input/output. This coordinate system and Kordulla' s coordinate 
system can be applied to any arbitrary wing configuration and cap 
be addressed to the entire wing region. J^From the above comparison, 
it is seen that the coordinate system used by Cebeci, et al. [8] is 
the most appropriate for the case of addre ssing an assu med wing con- 
figuration ^Thel BLAY uses this coordinate system.. 


Using the code indicated at the bottom of Figure 1, the non- 
cross flow coordinate system is expressed as 

>ndZ l -\ 2h,f, 2 vos0dccd2 - dy 2 (]_) 

and the compressible turbulent flow boundary layer equation is 
expressed as 

continuous equation: 

~(p uh 2 s\nO)+~(.flwh l sinO) 

°* dz (2) 

^ +JL — . 

. 1 Qy {.pv h i A*; s in 0 ) ’ 0 

x — momentum equation: 


4 



( 3 ) 


rnidti pw du , __ On 
“ ; — x~. - -- — 4/31;— -pu Aicotg 

/ii dx hz Dj- oy 


+P& ~A" 2 esc 0 + p u wK i2 ~ ~ 


esc 2 # dP 
A i dJ? 


cottfcsctfdP d . dll » 


z — momentum equation pu dw 


Energy equation: 


, pw aw — a w 2 

. + T" HT +Al '-q7 -f >w A* cote 

Ai a**? hi OX o y 

, , r , COtflcSC/? d7* 

+ PP A) CSC0 r /7W U»A 21 

A i d x 

esc 2 /? dp <9 / dw — :% 

- H — (/* pv w } 

hz dz ey x oy } 


pu qji pw on : — on _ d j n an 

A] ox + hz ez pv oy ~ ay\p r ay 

+p (i- jr)\ ( ur+w t +2uw C oso ) 


'V'l /'j 


(4) 


(5) 


Here /w^pv + iPv' ; p is density; p is pressure; u,w,v are the velocity 
components in the x,y,z directions; H is total enthalpy; y is the 
coefficient of viscosity; P r is the Prandtl number. The geometric 
parameters and K ^ are given by the following equations: 


l 


)_^Ail 


A ^wI^tfe (A2COse) '^"J 


(6a) 


1 Id/ x dA* 

As — — — — {— (Aicosfl) - -- — 
AiAssintf \dz Ox 


(6b) 


K — 1 ^ (tr u_ * dfl \ i n ( i r l ^ 30 \ 

A 13— ““7. % “-(Aj-f - — — )•+ COSfllAsf *- — ) 
“10 l Ai cat' v hzdz 1 


(6c) 


A ' j) ^ siuo f (a '- + a~I! )+ cos ^ a ' 4 i H )} 


(6d) 


The Reynolds' stress, ^ V ‘w’ and V//’ appearing in (3T through (5) 
is evaluated by means of the Cebeci type vortex viscous two-dimen- 
sional algebraic model 113] . 


Boundary conditions are adhesive adiabatic wall conditions at 
the wing surface: 
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y — 0 ; u — w = v = o , (dii/dy)^^]] - o 


(7) 


The value of the non-viscous flow on the outer fringe of the 

boundary layer is given by: 

y = d ; u = u t (x,z). ui.= w e (x, z). E-=E t ^ 

The initial conditions must be given by x = 0 (leading edge stag- 
nation line) z = 0 (wing root) ; or by z = (wing tip) . In the 

case of BLAY, the "stagnation equation" [3] is used at x = 0; the 
infinite sweep back wing equation [3] by which the z — integral term 
is set at zero for z = 0 or z = z tip* Thereby, the initial res- 
pective profiles are calculated. 


In the actual calculation, the Cebeci type [3] is converted 
into equations (2) through (5) and the difference scheme is employed 
in the boundary layer coordinate system. For the boundary layer 
conversion, not only can the peculiarity of the equation at x = 0 be 
eliminated, but in this type three-dimensional laminar flow/turbulent 
flow region, the boundary layer conversion is effective also from the 
standpoint of conserving the computer memory capacity and computer 
line . 


2 . 2 Numerical calculation method 

The three-dimensional boundary layer equation differs greatly 
in its character from the two-dimensional boundary layer equation 
in that the hyperbolic/parabolic forms, which are known as the basic 
principle of Raetz's zone of influence and zone of dependence, 
accrue to it. As shown in Figure 2, the disturbance generated at a 
point P within the boundary layer propagates into the normal line 
A-B on the body surface, having passed P instantaneously by virtue 
of its character as a parabolic type. Simultaneously, by virtue of 
its character as a hyperbolic type, it proceeds downstream along all 
stream lines passing through the line A-B. Therefore, the distur- 
bance generated at point P is transmitted to the entire wedge (cunei- 
form) region (the zone of influence) formed by the surface including 
the outermost streamline passing through the line A-B on the down- 
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Figure 2 . Three-dimensional boundary layer zones of 
influence and dependence. 

stream side. Conversely, the condition at point P is determined by 
the total reverse wedge region (zone of dependence) surrounded by 
the surface that includes the outermost stream line(s) passing 
through the line A-B on the upstream side. Therefore, when cal- 
culating the value of point P, it is impossible to perform a stable 
and reliable calculation unless a difference scheme is applied that 
will take into account all information on the zone of dependence. 

The numerical calculation method can be formulated by a variety 
of schemes by considering the character of all the above described 
boundary layer equations and depending on how to perform difference 
approximation in terms of differentials. The boundary layer equa- 
tions contain the following differential terms: 

<?( ) 5(] on d .di U (9) 

dz . oy . oy d y 

Using the symbols C3,0, +, and x to represent the above mentioned 
terms diagrammatically , Figure 3 shows the comparison between the /211 

various difference schemes and the BLAY shown in Table 1. As to the 
remark"^ on the Box method in Figure 3, it is to be noted that since 
the second order differential terms on y are converted to simultaneous 
first order differentials, the x symbol does not appear. Consider- 
ing in general the representative-difference schemes (a) through (e) 
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listed in Figure 3, it is seen that they are characterized by a 
number of properties as follows: (1) considering the stability of 

the difference scheme, it is implicit relative to the y direction; 
however, the interative process is introduced to linearize the non- 
linear terms and the computer operation load is increased, making 
the process inefficient [schemes (a) through (e) ] ; (2) the stability 

requirements on the numerical calculation dictate the need for a 
scheme that depends on the direction of flow. In other words, the 
difference configuration changes in response to the w code [schemes 
(a), (b) , (c) , (e) ] ; (3) being a totally gradual processing algorithm, 

a parallel processing operation cannot be performed [schemes (b) , (c) , 

(d) , (e)]; (4) the coordinate width in the z direction changes and 
when w is less than 0, the accuracy in the z direction deteriorates 
to the first order of coordinate width [schemes (a) through (e) ] . 


Property (2) described above is generally not suited to parallel 
processing operation because of the entry of the IF text or the pro- 
cedure corresponding to it within the computation loop. Also, in the 
presence of property (3) , the parallel processing operation is 
totally impossible. Calculation aerodynamics is a design tool and 
the computer to be used in the future will be a super computer for 
exclusive scientific and technical application which must have a 
parallel processing operation system as its basis. 


The BLAY difference scheme (predictor/corrector type Crank 
(algebraic) -Nicholson (PC-CN) scheme [11,13]) is a newly developed 
difference scheme which provides for higher efficiency and higher 
accuracy as well as being suitable for parallel processing operation 
by a specialized scientific and technical super computer. We next 
introduce this system briefly. 


The boundary layer equation considers the pressure gradient as 
functions of x, z whose characteristics are represented by a scalar 
equation that can be expressed as follows: 


dv. , 

a func[pc t z 

dx 


du. du \ d /,, 


(10) 
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Here, i,j,k represent the directional coordinates in the x,y,z 
directions, respectively; Ax., Ay., Az, represent the coordinate 

1 J K 

width immediately following those coordinates designators, and 
tabulate as u^k = u(x.,y.,z, ). The difference operators 6 and A 

d 13 k y y 

are formulated, respectively as follows: 




( 11 ) 


uLu-u'p 

J .v ( 4 >- u j* 5 =; 7y~ i+jyj ^ * jjij 


(12a) 




•O.-i" 


^±4*=y(^ + *s’ ± u> 


(12b) 


The difference operator 6^ attempts to become equation (11) . 
At this time, the PC-CN scheme is a two-stage semi-implicit scheme 
which can be expressed as follows (Figure 3, scheme (f ) ) : 
predicator : 


u t - - — +func(x :,•+ 1 , 2 * 

Jor,-/2 


«}*) 


(13) 



Corrector: 


— i+i 
i+| “>» 


4x, 


u vi , 1+4 1+4 

• /uiur( x I+ | ,z k ,u jk ", 8 X Ujk , 




2 


, / I + 1 \ y 1+7 a \ 

~4v u jk H A y ( J y U jk ) 


(14) 
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Figure 3. Difference schemes applied to three-dimensional 
boundary layer equation 

This scheme is characterized as follows: (1) It is linear 

with respect to unknown quantities (predictor step is u »4 s t ; 
corrector step is ti'* 1 ). All calculations are accomplished in 

two stages--predictor and corrector — with no iterative process. 

Therefore, as a whole, the computer operation load per coordinate 
point is low and the system is efficient; (2) since the z differ- 
ential term is expressed in a three-point difference core, the scheme 
itself depends on the w code and does not change; (3) as to the 
unknown quantity, it is implicit in the j direction, but explicit in 
the k .direction and mutually independent. Therefore, in the k 
direction array, parallel processing becomes possible; (4) the /212 

accuracy in the z direction is second order for Az; (5) evaluation 
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of the turbulent vortex viscosity can be performed directly at 1 + 

1/2 without iteration, thus facilitating the introduction of various 
high order turbulent models. Furthermore, it allows for consideration 
only of the CFL condition (Ax/Az less than u/w) on the x-y plane, as 
the stability condition for the scheme in the case of boundary layer 
equation in a manner similar to ther implicit schemes. 

2 . 3 Features of the BLAY 


Compared with other existing programs the BLAY possesses the 
following features: 

(1) the required computing time is the shortest 

(2) it is suitable for parallel processing computers 

(3) it has a second order accuracy in the z direction, too. 


As to item (1) above, there is no iterative convergence process 
in the respective steps; and the linearization of the ion linear 
terms can be traced to the fact that the process is simple and the 
operation load is low. We will next express this specifically with 
the example of the nonlinear term u.(ou/ox) . If, as a typical example 
of other high efficiency difference methods we cite the scheme that 
iteratively linearizes the second order accurate Newton Raphson 
model, it will be differentiated as follows (subscripts jk omitted) 


1 , l + l t\ 1 , i + l i) 

u d^ Y {u ^ ] YY.- (U ~ v ■ 


(~i+l , ~i + l i+l* / * \ ~ 

[u *(2 u --U )-(u . ) 


Here |7 ,+ 1 will be established as the initial estimated value 
or the pre-iteration value. We think of multiplication as a necessary 
operation load. If we assume the number of iterative cycles as the 
minimum 2, equation (15) becomes, for each coordinate point 

V 

multiplication load >2 x (iterative calculation: 

3 cycles) +1 = 7 (cycles) (16) 
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Even when we consider Cebeci ' s Box scheme, which is known as another 
high efficiency difference method and anticipating a minimum opera- 
tion cycle load by aggressively conserving an intermediate calcula- 
ted value, we obtain a value the same as equation (16). On the 
other hand, for the BLAY PC-CN scheme, as indicated by equations 
(3) and (14), they become 


(predictor) : 

du . 2 . ,+Y 

it—— — » u ■ (u *— 1 1 ) 

Ox jXj 

(17a) 


(corrector) : 

,, du. rs i , i +1 „ 
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"• *(tt -a } 

ox j.r 

(17b) 



Therefore, the number of multiplications for each coordinate point 
becomes : 

multiplication load =2+2=4 (cycles) (18) 

Let us now make a comparison using the actually measured value 
of the required computing time. To assure fairness, let us take the 
input/output portion of the BLAY without change and discuss the 
results of the numerical experiment conducted by replacing only the 
difference scheme. Taking as an example the problem of a tapered, 
backward swept wing configuration with coordinate grid numbers 1 x j 
x k x 2 (upper and lower wing surfaces) equal to 40 x 35 (average) 
x 41 x 2, and utilizing the FACOM M380 general purpose computer, the 
required computer time (excluding input/output) was 25 seconds for 
the linearized Newton Raphson model integrated into the Krause zig- 
zag implicit scheme; and 10 seconds for the PC-CN scheme. 

4 . Examples of application computation 

In order to conduct a computation, it is necessary as a bound- 
ary condition on the outer fringe of the boundary layer for the 
value b’f nonviscous flow to be given in terms of velocity components. 
The calculation results described below are derived, using non- 
viscous flow data on the outer fringe of the boundary layer as follows; 
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in the case of the low velocity incompressible flow problem, the 
results from the panel method; and for all other cases, the calcula- 
tion results from applying the transonic total potential flow anal- 
ysis program (AFPWING) [14], Here under the pressure gradient of 
the given boundary layer outer fringe the focus is on the feasibility 
of the BLAY being indicated as well as to denote what type of three- 
dimensional boundary layer flow is formed on the wing surface under 
that type of pressure distribution. Nonviscous-viscous interference 
is important, but is a separate problem, hence, will not be addressed 
here . 


Firstly, the results of computation in response to Brebner and 
Wyatt's experience [15] which was conducted for the purpose of 
investigating the validity of the BLAY boundary layer computation 
method are shown in Figure 4. In the figure, comparison is made only 
at one point but it is about the same even when several points are 
investigated. From this figure, it can also be seen for the boundary 
layer internal specific value that at the present time, even though 
many unsolved problems are yet encountered, it will be necessary to 
take into account viscous/nonviscous interference. Further, it will 
be necessary to investigate the turbulent flow model. 

Next, we show a calculation sample from the transonic flow 
problem. Generally, as a characteristic of transonic flow, in the 
reverse pressure gradient region, where the shock wave effect is 
strong even though localized, the boundary layer there peels off and /214 
the capability to calculate thereafter becomes ineffective. Since 
the BLAY scheme by virtue of its purpose prioritizes on the calcula- 
tion over the flow field in entirety, it has incorporated the FLARE 
[16] approximation which permits approximate calculation of small 
separation regions. This approximation is based on an estimate that 
the convection term (ud( )/Dx) within the reverse flow region is zero. 
Therefore, for a small, localized separation region it is known that 
a significant effort is not generated. However, in reality, there 
are many cases whe re computation can be performed without 
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Figure 4. Comparison of computed value and experimental 
value 
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Figure 5. Pressure distribution and boundary layer for 
the 720211 wing (wing upper surface). 
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Figure 6. Pressure distribution and boundary layer for the 
720211 wing (wing lower surface) 


incorporating a FLARE approximation by the fact that the shock wave 
itself is weakened numerically, and in the fact of applying turbulent 
flow of high Reynold's number. Results from four types of cases are 
shown below. Figures 5 through 9 show the results, graphically, of 
five sets of five plots which are from left to right: (a) pressure 

contours, (b) pressure distribution, (c) comparative plots of the 
body surface streamlines (limiting streamlines and broken lines) of 
the boundary layer outer fringe streamlines (continuous lines) and 
boundary layer, (d) boundary layer displacement thickness (A*) con- 
tours, (e) boundary layer displacement thickness distribution. The 
flow direction in all figures is from left to right. The terms in 
the figures are defined as follows: ALPHA: angle of attack; Re: /215 

Reynold's number; MACH: Mach number; YAW: yaw angle. Calculations 

were conducted for the four cases discussed below. We assume that 
the flow field is turbulent in practically all regions, and that 
laminar flow transitions to turbulent flow at x = 1% chord. 

•V' 


Firstly, the computation results on the transonic wing (Kamiya 
720211 wing [14]) are shown in Figure 5 (upper surface of wing) and 
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Figure 7. Pressure distribution and boundary layer over 
isolated swept back wings. 

Figure 6 (lower surface of wing) . The main stream Mach number is 
0.792; and the principal feature of the pressure distribution is 
that a strong shock wave exists on the upper wing surface's outer 
wing, and a pressure flattening exists in the vicinity of the wing 
root area. On the wing lower surface, the flow is subsonic over the 
entire region. The coefficient of lift, C_ , for this case is 0.566. 
We will first discuss the results of the boundary layer calculation 
along the wing upper surface. The calculation encounters the first 
reverse flow region (peeling off at the trailing edge) at 73% semi- 
span and 82% chord position, and the calculation continues thereafter 
to 86% chord by means of the FLARE approximation. However, since the 
separation region is large, the calculation thereafter becomes 
impossible. A special feature of the boundary layer that can be 
cited is that as seen in Figure 5 (d) a region of extremely strong 
three-dimensional characteristics appears in the wing root area. 

This type of region can also be observed in the calculation results 
of McLean's Boeing transport aircraft transonic wing. This region 
does not appear even in the case, for example, of calculation over 
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Figure 8. Pressure distribution and boundary layer over 
isolated forward swept wings. 
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Figure 9. Pressure distribution and boundary layer over 
isolated oblique wings. 
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a wing configured in a flat surface space distribution the same as /215 
for the 720211 wing with a NACA 0012 wing profile [17] . We will 
next discuss the lower surface of the wing (Figure 6) . The calcula- 
tion is performed to the trailing edge. The principal feature of 
the boundary layer is again that a region with a strong three-dimen- 
sional character is indicated in the wing root area. In the outer 
wing contours of displacement thickness appear parallel to the lead- 
ing edge, and two-dimensionally except at the trailing edge area. 

Next, let us consider three cases of independent wings (back- 
ward swept, forward swept and oblique) having a yaw angle in the 
flow field, and perform a comparative investigation of their boundary 
layers. Analyzed are planar configurations as shown in Figures 7 
through 9, whose cross sectional profiles are similar to the ONERA M6 . 

The three types of wing planar configurations are compared against a 
baseline configuration which provides a complete match in terms of 
the aspect ratio, the leading edge sweep back angle (±30°), and the 
chord length at various span stations. Conditions are: Mainstream 

Mach number, 0.92; angle of attack, 0°; yaw angle, 10°. First, let 
us compare the characteristics of the external nonviscous flow 
(Figures 7 through 9 (a) , (b) , (c) ) . The swept back wing acquires an 

additional yaw angle, resulting in a high, effective Mach number — 
in the case of the swept back wing (Figure 7) , on the left side; and 
in the case of the forward swept wing, on the right side. As a 
great difference in the respective pressure distributions, we can 
cite the existence of non-existence of a shock wave at the wing tip 
areas. The existence of a strong shock wave can be acknowledged for 
the swept back wing; but for the forward swept wing the status is 
"no shock wave" . For the oblique wing, because of the yaw angle, /217 
an additional 10° sweep back develops at the leading edge. As a 
result, the effective Mach number is low, and the condition over the 
entire wing area is "no shock wave" . The boundary layers under these 
pressure distributions are shown in Figures 7 through 9 (c) , (d) and 
(e) . First, considering the boundary layer surface streamlines, 
what is unique is that the flow (frictional resistance) for the swept 
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back wing is in the direction of the wing's external side and for 
the forward swept wing, in the direction of the wing's internal side. 
At the same time, the surface streamlines turn substantially at the 
point where the shock wave exists. On the other hand, in the case 
of the oblique wing, the surface streamlines in their entirely flow 
in the direction of the right wing tip which corresponds to the down- 
stream side. Next, comparing the boundary layer according to the 
displacement thickness in the case of the swept back wing, the dis- 
placement thickness is greatly increased at the shock wave position 
as can be seen in Figure 7 (d) and (e) , and the contour lines become 
more dense. 

On the other hand, if we consider the same aspect for the for- 
ward swept wing, it is seen from Figure 8(d) that the increase in 
the displacement thickness at the shock wave position is not as great 
as in the case of the swept back wing. If we consider the boundary 
layer in entirety, the displacement thickness distribution of the 
forward swept wing can be said, as shown in Figure 8 (d) to have its 
contours run in parallel with % chord in practically all regions, 
and two-dimensional. For the same aspect, in the case of the oblique 
wing, -it is even more two-dimensional (Figure 9(d)). Shown in Figure 
10 is a comparison of the distribution of displacement thickness at 
the central 50% span position (denoted by [> in Figures 7 through 9 (e) ) 
for the purpose of comparing the conditions surrounding the develop- 
ment of the boundary layer. Up to the 30% chord point, the boundary 
layer development is about the same for the swept back and forward 
swept wings; but beyond that and up to 55% chord, the displacement 
thickness increase is greater for the swept back wing. The rapid 
increase of the displacement thickness at the shock wave position is 
more rapid in the case of the swept back wing than with the forward 
swept wing proportionally to the strength of the shock wave. On the 
other _hand, the increase of the displacement thickness in the case 
of the oblique wing is monotonic as a whole, and merely increases 
abruptly in the trailing edge area. If we estimate the displacement 
thickness as a viscous effect, we see that the displacement thickness 
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increase is greatest for the oblique 
then the forward swept and swept 
back wings, in that order. 

5 .. Conclusions 

In this paper, we emphasized 
the BLAY's high efficiency of the 
boundary layer analysis program as 
an aerodynamic design tool for the 
viscous and nonviscous weak inter- 
ference problem (s) as a boundary 
layer correction in particular, to 
the transonic nonviscous flow analysis program, we find that it can 
be repeatedly used several tens of times. Therefore, its high effi- 
ciency is a most important element from cost considerations. 

The BLAY is a code developed from a baseline research code 
developed in 1980 by one of the authors whereby the input/output 
segment was generalized and a compensatory function added to the com- 
puting segment to achieve a pilot code, and then a production code. 

The BLAY is being applied to various flow problems, thus allowing for 
experience to be accumulated in its application scope, robustness, 
reliability, etc. What is introduced in this paper includes one 
portion (of that experience) , but it can be said that for the calcul- 
ation of boundary layers, they have in all cases been performed under 
rigorous conditions. It is believed that this program rates high in 
its reliability because of the simplicity of the scheme and its 
algorithms . 

We express our gratitude for the fact that Ms.Michiko Kikuchi 
of Fujitsu Incorporated (Kabushiki Kaisha) bore a large portion of 
the generalization work in the development of the BLAY. Mr. Kohei 
Suetsugu also of Fujitsu provided his good offices to expedite the 
use of the "AFPWING" program. We express our deepest thanks. 



Figure 10. Comparison of 
displacement thickness at 
the central areas of swept 
back, forward swept and 
oblique wings. 
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